import pandas as pd
import numpy as np
from math import pi
from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt
from scipy import signal
# from My_PCA import PCA as mypca
from sklearn.decomposition import PCA

if __name__ == "__main__":
    """
    从表格中将数据获取出来并做类型转换 
    """
    cols = range(10, 70)
    # print(cols)
    # dat = pd.read_csv("huxi 9.5 30.csv", usecols=cols)
    dat = pd.read_csv("fuwocheng 6.csv", usecols=cols)
    # dat = pd.read_csv("dunqi.csv", usecols=cols)
    # dat = pd.read_csv("jingzhi.csv", usecols=cols)
    data = np.array(dat)
    i = 0
    complex_data = []
    for m in data:
        complex_data.append([])
        for n in m:
            complex_data[i].append(complex(n))
        i += 1
    # print('complex_data[30][1] type is %s',type(complex_data[30][1]))# get the data type of complex_data list
    # print('"1+2j" type is %s',type(1+2j))
    # print(complex_data)

    """
    计算每一个数据的相位，若实部为0的相位为无穷大
    """
    i = 0
    phase_data = []
    phase_angle_data = []
    # print(type(data[0][0]))
    # rr=np.arctan(4/3)
    # mm=math.tan((60/180)*pi)
    # ii=2-float('inf') # float('inf')表示无穷大
    for m in data:
        phase_data.append([])
        phase_angle_data.append([])
        for n in m:
            real_one = complex(n).real
            imag_one = complex(n).imag
            if real_one == 0:
                phase_data[i].append(float('inf'))
                phase_angle_data[i].append(np.angle(complex(n)))
            else:
                phase_data[i].append(np.arctan(imag_one / real_one))
                phase_angle_data[i].append(np.angle(complex(n)))
        i += 1


    def sub_tow_antenna_data(indata):
        flag = 0
        sub_out_data = []
        for row_data in indata:
            sub_out_data.append([])
            for i in range(len(row_data) // 2):
                sub_results = row_data[30 + i] - row_data[i]  # rowdata
                sub_out_data[flag].append(sub_results)
            flag += 1
        return sub_out_data


    """
    将两组天线的相位角(arctan)也做减法存起来,存放到phase_sub_data中，并且存储方向跟表格方向一致，得到的数据个数是30组
    """
    phase_sub_data = sub_tow_antenna_data(phase_data)
    """
    将两组天线的相位角(angle)也做减法存起来,存放到phase_sub_data中，并且存储方向跟表格方向一致，得到的数据个数是30组
    """
    phase_angle_sub_data = sub_tow_antenna_data(phase_angle_data)
    """
    将两组天线发送过来的数据相减得到一组新的数据，两个天线各30个子载波，相减得到做减法后的数据集
    并且将相减后得到的30组数据的实部和虚部都提取出来
    
    """
    sub_data = []
    real = []
    imag = []
    r_flag = 0
    for rowdata in complex_data:
        sub_data.append([])
        real.append([])
        imag.append([])
        for i in range(len(rowdata) // 2):
            sub_results = complex(rowdata[30 + i] - rowdata[i])  # rowdata
            sub_data[r_flag].append(sub_results)
            real[r_flag].append(sub_results.real)
            imag[r_flag].append(sub_results.imag)
        r_flag += 1
    # print(complex_data)
    # print(sub_data)
    # print(real)
    # print(imag)
    """
    转置以后的数据,转置的目的在于从csv文件来的数据一个通道的数据是一列，转置后存储为一行
    是一个通道的数据，比如该段程序下方的phase_ready_data[0]就是第一通道的数据，
    phase_ready_data[0][0]就是第一通道第一个数据
    """
    real_size = np.shape(real)
    cnt_r = real_size[0]
    cnt_c = real_size[1]

    real_ready_data = []
    imag_ready_data = []
    complex_ready_data = []
    amp_ready_data = []
    phase_ready_data = []
    phase_ready_angle_data = []
    for c in range(cnt_c):
        real_ready = []
        imag_ready = []
        complex_ready = []
        amp_ready = []
        phase_ready = []
        angle_ready = []
        for r in range(cnt_r):
            real_temp = real[r][c]
            imag_temp = imag[r][c]
            amp_temp = np.sqrt(imag[r][c] ** 2 + real[r][c] ** 2)
            complex_temp = complex(real[r][c], imag[r][c])
            phase_temp = phase_sub_data[r][c]
            angle_temp = phase_angle_sub_data[r][c]

            real_ready.append(real_temp)
            imag_ready.append(imag_temp)
            amp_ready.append(amp_temp)
            complex_ready.append(complex_temp)
            phase_ready.append(phase_temp)
            angle_ready.append(angle_temp)

        real_ready_data.append(real_ready)
        imag_ready_data.append(imag_ready)
        amp_ready_data.append(amp_ready)
        complex_ready_data.append(complex_ready)
        phase_ready_data.append(phase_ready)
        phase_ready_angle_data.append(angle_ready)
    # print(real_ready_fft)
    # print(imag_ready_fft)

    """
    np.nan_to_num的作用是去除数据中的无效数据，
    """
    # print(real_ready_data[0])
    # print(type(real_ready_data[0]))

    real_array_data = np.nan_to_num(np.array(real_ready_data))
    complex_array_data = np.nan_to_num(np.array(complex_ready_data))
    amp_array_data = np.nan_to_num(np.array(amp_ready_data))
    phase_array_data = np.nan_to_num(np.array(phase_ready_data))
    phase_angle_array_data = np.nan_to_num(np.array(phase_ready_angle_data))
    """
    对数据进行滤波，用butter滤波器,滤波后的数据以filter标识  wn频率(Hz)/奈奎斯特频率（采样率*0.5）
    """
    fs = 100
    nyq = 0.5*fs
    lowcut = 1#1hz
    highcut = 1.4#1.5hz
    low = lowcut/nyq
    high = highcut/nyq
    [b, a] = signal.butter(3, 0.02, 'low')
    [b1, a1] = signal.butter(3, [low,high], 'band')
    amp_filter_data = signal.filtfilt(b, a, amp_array_data)
    real_filter_data = signal.filtfilt(b, a, real_array_data)
    complex_filter_data = signal.filtfilt(b, a, complex_array_data)
    phase_filter_data = np.nan_to_num(signal.filtfilt(b, a, phase_array_data))
    phase_filter_angle_data = np.nan_to_num((signal.filtfilt(b, a, phase_angle_array_data)))
    amp_bandfilter_data = signal.filtfilt(b1, a1, amp_array_data)
    """
    PCA进行主特征提取,貌似自己写的主特征提取不好使
    """
    pca = PCA(n_components=1)
    pca_data = pca.fit_transform(np.transpose(amp_filter_data))

    print('sklearn :\n{}'.format(pca_data))
    print(pca.explained_variance_ratio_)
    # pca_amp_data=mypca(amp_filter_data)
    # print("mypca:\n")
    # print(pca_amp_data)
    # size=np.shape(pca_amp_data[0])[1]
    plt.figure(1)
    # plt.subplot(311)
    # plt.plot([i for i in range(size)],np.transpose(pca_amp_data[0]))
    # plt.title("mypca")
    plt.subplot(211)
    #plt.plot(amp_filter_data[11])
    plt.plot(amp_bandfilter_data[11])
    plt.title("band channel 11")
    plt.subplot(212)
    plt.plot(pca_data)
    plt.title("sklearn  pca")
    # xx=range(len(pca_amp_data))
    # plt.plot(xx,-pca_amp_data)
    plt.figure(2)
    plt.plot(phase_filter_angle_data[11])
    plt.title("angle")
    """
    ！！！求方差最大并不能找到最合适的通道！！！还不知道什么方法好
    找amp_filter_data中的最高方差的通道
    #find the max secondmax thirdmax var
    #取三个最大值然后找到其所对应的索引（即通道）
    #max_index is the index of max varience channel
   
    """
    # amp_var = []
    # for amp_filter_data_var in amp_filter_data:
    #     amp_var.append(np.var(amp_filter_data_var))  # 求出所有通道方差并且存储到amp_var中
    #
    # sort_var = np.sort(amp_var)[26:30]
    # max_index = []
    # index = 0
    # for index in range(0, len(amp_var)):
    #     if amp_var[index] in sort_var:  # 找到方差等于最大四个值的索引（即第几个通道）
    #         max_index.append(index)
    # print(max_index)
    # print(amp_var)
    # # print(np.max(amp_var))
    # # print(np.where(amp_var == np.max(amp_var)))
    # var = np.var(amp_filter_data[max_index[0]])
    # var1 = np.var(amp_filter_data[max_index[1]])
    # var2 = np.var(amp_filter_data[max_index[2]])
    # var3 = np.var(amp_filter_data[max_index[3]])
    """
    求方差最大的四个通道的均值(即选出最好的四个通道求其均值作为分析使用)
    """
    # average_data = (amp_filter_data[max_index[0]] + amp_filter_data[max_index[1]] + amp_filter_data[max_index[2]] +
    #                 amp_filter_data[max_index[3]]) / 4
    # # average_data =  amp_filter_data[max_index[2]]
    # # 极大值点
    # print(average_data[signal.argrelextrema(average_data, np.greater)])
    # print(signal.argrelextrema(average_data, np.greater))
    # x_maximumpoint=np.transpose(np.array(signal.argrelextrema(average_data, np.greater),dtype=int))
    # # 极小值点
    # print(average_data[signal.argrelextrema(-average_data,np.greater)])
    # print(signal.argrelextrema(-average_data,np.greater))
    # print(type(signal.argrelextrema(average_data,np.greater)))
    # x_minimumpoint = signal.argrelextrema(average_data, np.greater)
    """
    画图部分
    """

    plt.show()
